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Abstract. We develop a analytical and quantitative theory of the tube model concept for entangled net- 
works of semiflexible polymers. The absolute value of the tube diameter L± is derived as a function of the 

polymers' persistence length l p and mesh size £ of the network. To leading order we find L± = 0.32£ 6 / 5 Zp 1//s , 
which is consistent with known asymptotic scaling laws. Additionally, our theory provides corrections to 
scaling that account for finite polymer length effects and are dominated by the mesh size to polymer 
length ratio. We support our analytical studies by extensive computer simulations. These allow to verify 
assumptions essential to our theoretical description and provide an excellent agreement with the analyt- 
ically calculated tube diameter. Furthermore, we present simulation data for the distribution function of 
tube widths in the network. 

PACS. 83.10.Kn Reptation and tube theories - 87.15.Aa Theory and modeling; computer simulation - 
87.16.Ka Filaments, microtubules, their networks, and supramolecular assemblies 



1 Introduction 

Filamentous actin (F-actin) is a semiflexible biopolymer 
that has been the object of intensive research from sev- 
eral domains. As a major constituent of the cytoskeleton, 
F-actin networks play a key role in the ongoing puzzle of 
cell mechanics [TJ[2] and cell motility 3 . Depending on the 
presence of binding proteins, F-actin strands at medium 
concentration can form both chemical (cross-linked) and 
physical (entangled) networks with different elastic prop- 
erties [HOIS]. Besides rheological methods, single polymer 
visualizations also are feasible, as the strands exceed most 
synthetic polymers by length. This facilitated the observa- 
tion of tube-like regions along which DNA or F-actin fila- 
ments reptate [7][8] . The confinement of polymers to these 
cylindrical cages confirmed the "tube model" postulated 
earlier by de Gennes [9] and Doi and Edwards [10]. This 
long standing paradigm had proven a successful concept 
to reduce the complex structure of entangled networks to 
a single polymer problem. 

In such entangled networks, polymers can effortlessly 
slide past each other but are not allowed to cross. Their 
interaction is thus mainly of entropic nature as entangle- 
ments mutually restrict the accessible configuration space. 
Grasping this feature in a single polymer model has lead to 
the famous tube model [9lll0j. The suppression of trans- 
verse undulations of a test polymer by the surrounding 
polymers (Fig[T]) is modelled by a tube. This tube follows 
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Fig. 1. The effect of all surrounding polymers that hinder 
the test polymer's transverse displacement is described by a 
hypothetical tube. 



the average path of the test polymer and its profile is fre- 
quently modelled by a harmonic potential. The average 
strength of this potential is determined by the local den- 
sity of the network. The tube concept has proven a suc- 
cessful tool to derive scaling laws for several network prop- 
erties [TOirrT] . For example, due to the confinement energy 
of the filament inside, the tube diameter can be connected 
to mechanical properties of the network, e.g. the different 
moduli 6,12,13 . However, due to the phenomenological 
nature of the tube model, most of its benefits have been 
mainly qualitative. Recently, also quantitative predictions 
of the plateau modulus and the tube diameter of flexible 
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polymers melts were achieved by a novel approach based 
on the microscopic foundations and the topological struc- 
ture of the network [T4lfT5j . 

Even if most concepts developed for flexible polymers 
can not be carried over to the semifiexible case with its 
large persistence length, the tube model is perfectly appli- 
cable as well. However, while in general scaling laws of the 
tube diameter [16] or the plateau modulus [13] are well es- 
tablished, quantitative theories are still under debate and 
lack from approval by measurements of sufficient accuracy. 
Again the challenge is to make the successful tube model 
quantitative by connecting the phenomenological tube and 
its microscopic origins. In the present work, we contribute 
to the discussion by supplying an absolute value for the 
tube diameter from a theory supported by extensive com- 
puter simulations. 

We will proceed as follows: in Section [2] the model un- 
der investigation is defined and all relevant length scales 
are discussed. By the analysis of the free energy cost for 
confining a polymer to a hypothetical tube, the tube diam- 
eter is derived as a function of Odijk's deflection length for 
finite-length polymers. The appropriate deflection length 
for a given polymer concentration and persistence length is 
derived in the following sections. To this end, the polymer 
is modeled by a sequence of independent rods in Section 
[H Criteria for the correct choice of the independent rod 
length and a self-consistent determination of the tube di- 
ameter are developed, before a final result for the tube 
diameter is obtained in Section 01 Extensive numerical 
simulations supporting these result and providing addi- 
tional insight are presented in Section [5] followed by our 
conclusion in Section [6j 



2 Model Definition 



We consider a mono-disperse network of physically entan- 
gled polymers with a particular focus on pure solutions of 
the biopolymer F-actin. The polymer density is given by 
the number v of polymers of length L per unit volume. 
The polymers are of bending stiffness k corresponding to 
a persistence length l p = n/hsT. A single polymer's con- 
figuration r(s) is parameterized by the arc length s and 
the average distance between the polymer chains can be 
characterized by a mesh size £ := y/3/(uL) B We will 
describe the constituent polymers by the worm-like chain 
model 17,18 and exploit the tube model concept [9 , 1 OJ 
to reduce the description of the network to a single poly- 
mer and its neighbors. In the following we will begin our 
analysis with an investigation of the different length scales 
involved in the system. 



1 The mesh size has the unit length and can be interpreted as 
an average distance between network constituents. While the 
denominator yJlfvL ensures the correct scaling the numerator 
is a mere definition. 



2.1 Length Scales 

Typical F-action solutions are polydisperse with a mean 
filament length L « 22/xm [19] , With a persistence length 
l p » 17//m [2P1I2T] comparable to its length, it is the text- 
book example of a semifiexible polymer. At a concentra- 
tion of c = 0.5 mg/ml corresponding to v sa 1/im -3 |22j 
the average mesh size equals £ ~ OA/im. We can thus 
state that the persistence length of a filament is much 
larger than the distance to its neighbors, l p ^> £. Since 
the tube diameter L± is at most of the order of the mesh 
size, this additionally implies l p Lj_. The polymer will 
thus not deviate far from the tube center. Consequently, 
configurations where the polymer folds back onto itself are 
rendered unlikely. This is a minimal requirement to model 
the tube by a harmonic potential of strength 7. The po- 
tential has to be seen as a hypothetical tube representing 
the joint contribution of all surrounding polymers which 
constrain the transverse undulations of a given polymer 
(see Fig H]}. 

The energy of a certain polymer contour r(s) is the 
sum of the bending energy of the polymer and its con- 
finement into the harmonic potential and is given in the 
weakly-bending rod approximation by 



H( 7> k) 



ds 



(1) 



Here r(s) = (s, r±_ (s)) is a parameterization in arc-length s 
and transverse displacement rj_(s) = {y(s), z(s)) from the 
tube center. The prime denotes a derivative with respect 
to s. This harmonic approximation to the Hamiltonian 
of the worm-like chain model is valid as long as |r" | <C 
1, i.e. as long as the transverse coordinates of the tube 
coordinate can be considered to remain single valued. 

With the thermal average (•) the tube diameter can 
now be defined as 



L 



1 



( / dsr 2 ± (s)) 



(2) 



So far we have identified two length scales: the length scale 
of persistence length and the total polymer length describ- 
ing the properties of one specific polymer, and the length 
scale of mesh size and the tube diameter describing the 
properties of the network structure. Additionally we in- 
troduce the deflection length := as a third 
useful length scale. It is interpreted below as that length 
on which interactions between single polymer and network 
occur. More precisely, it is a measure for the number of 
contacts of the polymer with the tube walls. For large 
confinement strength 7 the tube is small, making interac- 
tion with the encaged polymer more likely and therefore 
resulting in a small deflection length. On the other hand, 
for a large polymer rigidity k transverse undulations allow- 
ing contacts with the tube walls are energetically unfavor- 
able and the distance between contact will decrease. For 
l p 3> L± we expect the deflection length to be distinctively 
smaller than the polymer length, but also larger than the 
tube diameter. For quantification we consider the free en- 
ergy cost AF(-f) of confining the polymer to the tube. It 
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can be found from the partition sum that is obtained as a 
path integral over all polymer configurations: 

exp [-/3AF( 7 )] = J T>[r ± (a)] exp[-/3H (k, 7)] (3) 

with [3 — l/k B T. In the limit of infinitely long polymers 
the free energy cost is [23 



AF = V2k B T 



L 



(4) 



This result fits into the picture of the deflection length as 
measure for the average distance between successive colli- 
sions of the polymer with its tube. If the typical distance 
between two collisions is given by L<j, the free energy loss 
results as the sum over all L/L^ points of contact where 
every collision costs one k B T. The free energy now allows 
one to derive the tube diameter as 



L 2 , = 



2 dAF 



T 3 



L d 1 V2/p 



(5) 



In the limit of infinite polymer length we have thus derived 
the tube diameter as a function of the deflection length by 
differentiation of the free energy cost. 

The above consideration also sets the road map for the 
remaining work. To calculate the tube diameter for the 
network, we need first to connect free energy and tube 
diameter for polymers of finite length and then derive the 
deflection length for the model under investigation. 



2.2 Finite length Polymers 

For finite size polymers the path integral in Eq.([3]) can 
be evaluated exactly 23 ; 24J and with the dimensionless 
deflection length l& := L^/L results in 



AF = -2k B Tg(l d ) 



with 



. 9 (/ d ) = ln(^)-iln(sinh 2 7 l--sin 2 7 L- 



(6) 



(7) 



The limit of small Id that is guaranteed by L ^S> La as 
stated above, allows an expansion 

9 (Z d ) = - * +hi(Z2) + (9( e -V'd ) for ?d ^ , (8) 

where the first term is just the result for polymers with 
infinite length ([4]). Upon again using the relation L\ = 
(2/L)(dAF(l d )/dj) with the inner derivative dl^/d"/ = 
— (X 4 /4k)^ the tube diameter becomes 



Li 



(9) 



For later convenience we simplify this to l\ — h(l&) by 
introducing a dimensionless tube width lj_ and function 
h(x) as 



L ± 1 P 
L 3 



and 



h(x) :-- 



x 5 g'(x) 



(10) 



This relation connects the wanted tube diameter to the 
deflection length and hence to the hypothetical tube po- 
tential 7 at a given bending rigidity. In the following we 
will further investigate the tube properties and set up a 
model that allows one to derive the deflection length and 
thereby the hypothetical harmonic tube potential strength 
from the polymer concentration and persistence length. 



3 Independent Rod Model 

For simplification and as an anticipation towards the com- 
puter simulations, consider for the time being a polymer in 
a two-dimensional (2D) plane. In this case the transverse 
displacement vector r_L reduces to a single component. 
The undulations of the test polymer in 2D are hindered by 
point-like obstacles as depicted in Fig. [2] (top). These ob- 
stacles represent the cuts of the surrounding polymers in 
three dimensions with the chosen fluctuation plane. Given 
an appropriate number of 2D obstacles equivalent to the 
density of surrounding polymers in 3D, the transverse dis- 
placement will correspond to one of the two components 
of the displacement vector rj_, if we assume the fluctua- 
tions of these components to be independent. Bearing in 
mind the large persistence length compared to the mesh 
size, the surrounding polymers in 3D are modelled as rigid 
rods and the area density pmc of obstacles in 2D corre- 
sponding to a polymer concentration v in the 3D network 
is Pmc = 2vL/ir. It is computed in Appendix lAl and will 
be explicitly needed for the comparison with simulation 
results. 

Recalling the Hamiltonian JT]), the polymer's free en- 
ergy has a bending and an entropic contribution. To mini- 
mize the free energy it can be favorable to trade in bending 
energy for a wider tube. Thereby entropy is gained due to 
a larger available free volume, but the polymer is forced to 
sacrifice energy to obtain its curvature (see Fig. [5] (top)). 
This competition defines a characteristic length L that has 
to be of the order of the deflection length La, since this 
is the length scale characterizing interaction of the test 
polymer and its environment. 

In the following we will develop an analytical theory 
based on an independent rod model (IRM) that is in- 
spired by the competition we have just discussed. To this 
end, we use a simplified model of a semiflexible polymer, 
in which the flexibility is localized to the joints of a se- 
quence of independent stiff rods of length L. After deriv- 
ing the transverse fluctuations of a single independent rod 
in an environment of fluctuating neighbors, we apply a 
self-consistency argument to arrive at the corresponding 
tube width of the full length semiflexible polymer. Note 
that the analysis is carried out for three dimensions and 
the 2D simplification only serves for illustration and for 
simulations later on. 

To begin with, consider the test polymer to be divided 
into independent segments of length L that are assumed to 
be completely rigid rods and are only allowed to undergo 
transverse fluctuations. As the flexibility in the IRM de- 
pends on the number of joints, it is obvious that the choice 
of L is crucial for the resulting tube diameter. Picturing 
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Fig. 2. (top) The fluctuation tube of a semifiexible polymer in 
a network of constraints is determined by a delicate balance of 
cntropic and bending energy, (middle) Scheme of decomposi- 
tion of a semifiexible polymer into rigid rods of length L. The 
flexibility is localized to the joints between independent rods. 
Given the proper choice of L both models produce the same 
transverse fluctuation area, (bottom) Small rod length overes- 
timates and large rod length underestimates fluctuations. 



the decomposition of the test polymer as in Fig. [2] (mid- 
dle) it can be seen that the transverse fluctuations of the 
independent rods are hindered by the two closest obstacle 
polymers normal to either side of each segment of length 
L. If L is chosen too large (e.g. L = L in the worst case) 
the area of transverse fluctuations will be much smaller 
than for a true semifiexible polymer because flexibility is 
underestimated (Fig. [2] (bottom, right)). On the contrary, 
if L is chosen too small, the normal distance to the near- 
est obstacle can be quite large (Fig. [5] (bottom, left)). This 
overestimation of flexibility results in a transverse fluctu- 
ation area that is large compared to the polymer we try 
to model. Before we further discuss the proper choice of 
L, we will focus on the behavior of a single independent 
rod in more detail. 



The transverse fluctuation of a single stiff rod in the 
(y, z)-plane are constrained by the projections of the sur- 
rounding network constituents to this plane as depicted 
in Fig. [3] (left). Since the mesh size is much smaller than 
the persistence length, the surrounding polymers can be 
assumed to be straight and "dangling ends" are neglected. 
The size of the shaded cross section will decrease with in- 
creasing density of polymers, i.e. with a decreased mesh 
size. Thus the tube diameter is of order of the mesh size 
and scales as L± oc £ for a given 2D plane. Furthermore, 
an increase of the length L of the rigid rod signifies an in- 
crease of obstacles that will be projected to the plane. As 
the average distance between surrounding polymers in di- 
rection of the test rod is also given by the mesh size £, the 
average number projected onto the plane increases as £/£. 
As this reduces the cross section area, we finally arrive at 
an overall scaling of the tube diameter as L± cx £ 2 /L. 




Fig. 3. (left) Projection of constraining polymers to the plane 
of transverse fluctuations of a test polymer (black dot). As the 
mesh size is much smaller than the persistence length, the con- 
straining filaments, can be assumed to be straight. The shaded 
area is the accessible tube area for a specific obstacle configu- 
ration, (right) Corresponding setup for a simplified geometry 
where obstacles can only be aligned with coordinate axes. 

Before we quantify this scaling result in the next sec- 
tion, let us first have a closer look at the obstacles. In a 
self-consistent treatment these evidently have themselves 
to be regarded as semifiexible polymers of the network and 
therefore undergo fluctuations around an average position 
as well. This causes the cross section area to smear out, 
as the test polymer has now a non-vanishing probability 
to take on values behind the average obstacle position. 
In terms of a confinement potential the cross section is no 
more described by an infinite well, but by some continuous 
potential which earlier has been assumed to be harmonic 
with strength 7 per unit length of a polymer. The ob- 
stacle fluctuations will also be modelled as Gaussian and 
to distinguish between the test polymer mean square dis- 
placement L\ and the obstacle's, the latter is denoted as 
a 2 . In a self-consistent treatment of the network the aver- 
age tube width L± of the test polymer is then determined 
as a function of the obstacle fluctuations a, where a is cho- 
sen such that L± = a. Of course, the value L± of a single 
obstacle configuration will not only depend on a but also 
on the obstacle positions in that specific configuration. 
Consequently, averaging over all obstacle configurations 
will result in a distribution P(L±) and self-consistency 
would then also require a distribution P(cr). However, if 
we assume these distributions to be reasonably peaked, 
we can use their averages as a good approximation. The 
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self-consistency of distributed tube widths is verified by 
simulations in Section [5] 

3.1 Single stiff rod in simplified geometry 

According to the assumptions made above the obstacles 
(in a top view) are completely described by a normal dis- 
tance r k from the test polymer and an orientation a k ; 
compare Fig. [3] We will neglect correlations and assume 
the obstacles to be uniformly distributed. The probability 
to find an obstacle with a certain direction at a speci- 
fied point is independent of the direction and that point. 
This corresponds to a complete factorization of the net- 
work distribution function into single polymer distribution 
functions. 

Consider first a simplified geometry in which all obsta- 
cles are either parallel to the y or the z axis as depicted 
in Fig. [3] (right) . As fluctuations in both coordinates are 
assumed to be independent and equivalent, the task of 
computing the tube width is reduced to a one dimen- 
sional problem with a single coordinate r. The network 
density or mesh size enters as the number p of obstacles 
per unit length. This density should be chosen such, that 
the average number of obstacles at a certain distance r 
from the test rod in the IRM is the same as the average 
number of obstacle polymers featuring a minimal distance 
r from the test polymer. This density is proportional to 
the length L of the stiff segment and the number of sur- 
rounding polymers in a unit volume vL. The exact relation 
p = (n/2)(vLL/A) is calculated in Appendix [B] 

As the obstacles are assumed to undergo Gaussian fluc- 
tuations around their average position r k , the correspond- 
ing probability density is 

P (r-r k ,a) := (2™ 2 )~ 1/2 e^?' . (n) 

If the test rod interacts with only a single obstacle, we 
can state that the probability to find the test rod at a 
certain position is given by the fraction of realizations still 
accessible to the obstacle. In this case 

poo 

P+{r,r k ,a) = dr'P (r' - r k ,a) (12) 

J r 

is the fraction of configuration space still accessible to the 
obstacle if the test rod is placed at r (for r k > 0). Com- 
pleting the integral yields 

P+(r,r k ,a) = i e rfc (13) 

and the corresponding probability for obstacles at negative 
positions P_ (r, r k , a) is simply obtained by a inverted sign 
of the argument. The probability to find the test rod at 
a position r for a given configuration of obstacles {r k } is 
then given by the product of all probabilities 

P(r,{r k },a) = ^ J] P+(r,r k ,a) J] P -(r, r k , a) . 

k,r k >0 k.r k <0 

(14) 



The normalization N = N ({r k }, cr) is determined by the 
condition J drP(r, {r k } 7 a) = 1 and depends on the obsta- 
cle configuration. 

As the function P + (r, r k , cr) reduces to a Heaviside 
function in the case of cr — > 0, the product in Eq. (TT4")) 
can be written as 8(r — r_)#(r — r +)/( r + ~ r -) where 
r + and r_ are the positions of the two closest obstacles. 
This reduction is justified because all further obstacles 
are completely shadowed by the two nearest neighbors. 
In the case of a non-vanishing a the probability distribu- 
tion P(r, {r kl a k }, cr) will not be rectangular anymore but 
smear out. The test rod has a non- vanishing probability 
to be found behind the average position of the closest ob- 
stacle and thus a chance to feel the interaction of further 
network constituents. However, sketching the distribution 
in Fig. IU it becomes intuitively clear that this probability 
rapidly approaches zero for far obstacles or small fluctua- 
tion amplitudes a . We will exploit this fact in the numer- 
ical analysis below and in the simulations. 




-5-2 2 3 7 

transverse coordinate 



Fig. 4. Probability density to find the test rod at a spatial po- 
sition for mutual interaction with a single obstacle (solid lines) 
and resulting probability in an environment of all obstacles 
(dashed line). The x-axis tics mark the center position of each 
obstacle. Distant obstacle only have a negligible influence on 
the overall probability function. 

The distribution function, Eq. p^]) , for the test rod at 
hand, averages of any function f(r) can now be calculated 
for a single realization of obstacles as 



/( r W = J drf(r)P(r,{r k },<r) , (15) 

where the index {r k } denotes the specific obstacle config- 
uration. The tube center of the test rod is then 

r{{r k },a) ~r {rh} (16) 

and the width of the probability distribution is the wanted 
tube diameter 

L 2 ± ({r k },a) ;=^ {rk} -r\ rk} . (17) 
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The derived tube diameter of the test rod is not only a 
function of the fluctuation width a but also of the spe- 
cific obstacle configuration. Consequently, sampling over 
different obstacle sets will result in a distribution of values 
for L±. As mentioned earlier this distribution should be 
described by a single characteristic value - consistent with 
the obstacle fluctuations that have also be assumed to be 
of equal size. Since the obstacles are uniformly distributed, 
they can be fully described by the density encoded in the 
average number of obstacles p per line. This is achieved 
by integrating out all obstacle positions and orientations 
in L^({rfc},cr) to arrive at the only density dependent 
L\(p : a). We choose a simple average over a large number 
N of obstacle sets {r k } like 

/ N rR/2 , \ 

(/(k})>p = ( n j m j /«»•*}) . ( is ) 

where R — N/p. In this nomenclature the average tube 
diameter0is obtained as L\(p,a) = (L^({r k }, a)) p . 
Self-consistency is now expressed as 

L\{p,a) = a 2 (19) 

at the point of self-consistency (PSC) a = a* . By measur- 
ing length in 1/p we can rewrite this to a dimensionless 
master curve I (per): 

Li(p,a) = ±l(pa), (20) 
P 

since L±,l/p and a are all lengths, The task of finding 
the self-consistent tube width L±(p,a*) = a* translates 
to finding 1(C) = C 2 where the constant C — pa*. As 
soon as this is achieved, the self-consistent tube diameter 
is available as a function of density p only and hence it 
depends like 

1, = C - = ft (21) 

on the rod length L and mesh size £. 

The numerical determination of C is achieved by an in- 
tegration using a Monte-Carlo procedure. It includes the 
TV-fold integrals over the obstacle positions r k from Eq. 
(fT4"l) as well as the integration over the test polymer posi- 
tion r from Eq. (fT5|) . As mentioned above, the probability 
distribution rapidly decreases at distances far from the 
closest obstacle. Hence, we have restricted the integration 
range of the dr integration in the Monte-Carlo samples to 
values [y m in ~ 5ct, y max + 5cr] where y mm and y max are 
the closest obstacle at either side. Furthermore, the fast 
decrease of the probability distribution renders the contri- 
bution of distant obstacles quasi to zero. We can therefore 
drop all obstacles with y k ^ [y m in — 10c, y ma x + lOc]. The 

2 Of course, one could also image a different characterization 
of the average tube diameter, e.g. the median or the maximal 
diameter. We choose the average as the most obvious quantity 
experimental groups might measure, e.g. in analyzing different 
fluorescent microscopy images. 



results are depicted in Fig. [5] and a graphical solution for 
the PSC constant results in 

C « 3.64 . (22) 

Special attention should be paid to the behavior at pa = 0. 
It provides a good test whether the used IRM is adequate 
and allows for a verification of the numerics. At a finite 
density as required by the tube concept, 1(0) reflects the 
situation of immobile obstacles with a = 0. At this point 
the tube diameter should remain finite and its value should 
be given by the density of obstacles. From the obstacle 
statistics and density per unit length p/4, the probability 
to find the first obstacle at position r± is known to be 
P(r±) = exp(— r±p/i). In the case of fixed obstacles the 
available fluctuation area is 2Lj_ = r + —r- and the expec- 
tation value 4(L 2 L ) — ((r + — r_) 2 ) can be computed from 
the probability density above. Taking care of the normal- 
ization one arrives at L± = \/8/p. The master function 
yields 1(0) — p 2 L\ = 8, a value in good agreement with 
the data (circles) in Fig. [5l 

3.2 Generic 2d Geometry 

If the simplification of axis-parallel obstacle polymers is 
dropped again, the obstacle configuration needs to be spec- 
ified by a set of radii {rk} and angles The probability 
to find the test rod at a position (y, z) for a given config- 
uration of obstacles {rk,ctk} in the two-dimensional case 
as in Fig. [3] (left) is then again given by the product of all 
probabilities where different angles have to be accounted 
for: 

P(y,z,{r k ,a k },a) = — J| P± (y cos a k + z sin a k , r k , a) . 

k 

The normalization factor N = N({r k , a k }, a) is again de- 
termined by the condition J dydzP(y, z, {r kl a k }, a) = 1. 

In a single obstacle configuration the tube diameters 
L± y ,z in the y and z direction will in general be differ- 
ent. However, in averaging over all configurations isotropy 
must be recovered to show 

L\(p,a) = ( J Li y ({r /c ,a fc },cr)) p = (L\ z ({r h , a k }, a)) p . 

The average over obstacle configurations at fixed density 
of uniformly distributed obstacles is performed as 

(/({rw)> P = (n/ fl ir jT^t) fi{rk ' ak}) 

(25) 

with integration range being again R = N/p. Note that 
contrary to the simplified geometry the obstacle density 
per unit length p in this case is given by p — (n /2)(vLL). 
Evaluating the integrals in Eq. (|25[) again by the Monte- 
Carlo method results in the data plotted in Fig. [5] (tri- 
angles) , where suppression of irrelevant obstacles was im- 
plied analog to the simplified geometry. The results do not 
deviate much from the data obtained earlier (circles), i.e. 
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Fig. 5. Master curve l(pcr) of the tube diameter rescaled 
by obstacle density obtained by MC simulation for simplified 
(circles) and generic geometry (triangles); intersection with 
quadratic obstacle fluctuation amplitude marks the point of 
self-consistence. The error of the simplified geometry is sur- 
prisingly small. 



the mistake in using a simplified geometry is surprisingly 
small. Again the value of the PSC is obtained graphically. 
It yields 

C w 3.52 (26) 
and will be used in the remainder of this work. 



3.3 Choice of Independent Rod Length 

As discussed before and illustrated in Fig. [2] (bottom) the 
choice of L is crucial for the success of the IRM. The 
number L/L of independent rods can be regarded as a 
measure for the flexibility of the modeled polymer and 
has to be chosen such that the transverse excursions of 
the ensemble of stiff rods equal the fluctuations of the 
actual semiflexible polymer. To this end we consider both 
systems in a generic harmonic potential 



U[y(s)] 



7 



(27) 



with the potential minimum y (s) as a Gaussian variable 
with (y°(s)y°(s')} — aS(s — s'). This corresponds to the 
assumption of a "Gaussian random backbone" as a general 
property of the tube. We use this intuitive assumption as 
one possible prerequisite to determine the segment length 
L. Of course, other possibilities can be imagined. Note 
that the simulations in Sec. [5] will justify this assumptions 
a posteriori. 

The average position y(s) as a functional of a given 
potential y°(s) is obtained as an average over all poly- 
mer configurations in this potential. Averaging then over 
all potential conformations yields the the mean square of 

2 

the polymer's transverse fluctuations (y(s) ). The over- 
line thus denotes an average in a given potential and the 
brackets denote an average over all potentials. While the 



transverse fluctuations of a rigid rod are a function of the 
potential parameters a, 7 only the response of a semiflex- 
ible polymer will additionally depend on its stiffness. This 
evidently provides a tool to connect the semiflexible poly- 
mer persistence length and the length L from the IRM by 

2 

demanding that the fluctuations (y(s) ) for given poten- 
tial parameters a, 7 are the same for both cases. 

Starting with the IRM, it is sufficient to consider only 
one stiff rod, as the individual rods are statistically inde- 
pendent. The average position is then 



(28) 



and the transverse fluctuations 



a 



L d S £ ds >(y°( S )yV( S >)) = ^. (29) 



For the semiflexible polymer the fluctuations of polymer 
and tube potential are decomposed into modes (Appendix 
0: 



(30) 



where the mode analysis yields yu = 2^/(1 + it^d) ( com 



pare Eq. (|5"Tj) ) with q% w 7r(fc + 1/2). Using now the corre- 
lations of the Gaussian random tube profile and the iden- 
tity (f53|) the polymer fluctuations can be related to the 
deflection length as: 



(t) = 



1 



a h'(l d ) 
L 411 



(31) 



Equating the fluctuations for the IRM and the semiflexible 
polymer fixes the segment length to 



L = L 



h'(l d ) 



(32) 



Concluding the last section, we have obtained the tube 
diameter for a sequence of independent rods of length L 
and derived a condition how to fix this length to correctly 
mimic the behavior of a semiflexible polymer in a network 
of same mesh size. It has turned out the the criteria for 
the correct rod size is a function of the deflection length. 



4 Results 

If we recall that the tube diameter for a semiflexible poly- 
mer was derived in Sec. [2] from the Hamiltonian with 
a likewise dependence on deflection length, we are now 
equipped to set up an implicit equation to determine this 
deflection length. Afterwards the tube diameter can be 
derived from a simple calculation. 

Equating the expressions for the tube diameter of the 
polymer (JTDJ) and the IRM (|2"Tj) respectively yields 



L 



4C 2 ^ 



(33) 



8 



Hauke Hinsch et al.: Quantitative Tube Model for Semifiexible Polymer Solutions 



With the correct rod length (|32f the implicit equation for 
the dimensionless deflection length is 



h(l d ) 



c 2 [h'(i d )] 2 i P e 

32tt 2 ?« L 5 



(34) 



Solving this equation, determines l d from the system's pa- 
rameter l p , L and £ . It is achieved by introducing a dimen- 
sionless function 



Inversion then yields 



(35) 



(36) 



With the abbreviation D = 3C/ir we finally obtain for 
the dimensional deflection length to first order and second 
order in the argument of j~ : 

r>2/5 n4/5 t8/5/2/5 

T - ^4/5,1/5 | U ? 'P /o 7 <, 

^ d - 2 13 / 10 p 2 n /i03 L ' [ ' 
By application of (fTT))) . the tube diameter is easily obtained 

D 3/5 £6/5 D £2 



as 



2 27/10 ,1/5 2 5 / 2 L 
lp 



(38) 



Evaluation of the numerical factors holds the following 
results to first order 



L d « 0.66C 4/5 ^ /5 , £j 



0.32 



£6/5 

'p 



(39) 



Note that this also determines the confinement free en- 
ergy of the polymer to AF w 2.14 fc B T L 1/5 . The lead- 

ing term of the tube diameter agrees with the established 
scaling [16J . The additional term's dependence on the in- 
verse polymer length indicates a finite length effect. It 
can be traced back to the partition sum of a finite poly- 
mer © and accounts for boundary effects at the end of 
the tube. If the free energy of infinite polymers (fj| is used 
throughout the calculations, all higher order terms van- 
ish accordingly. In an earlier work [25 another prefactor 

of L± w 0.53£ 6 / 5 Z p l ' 5 for the scaling term has been pre- 
dicted by rather different accounting of obstacles. 

It is important to be aware of the subtle difference 
between the explicit length dependence of the first order 
term and the implicit dependence on L that enters via 
the mesh size £ = \fiJvL. In a polydisperse polymer so- 
lution the L in the mesh size has to be the average poly- 
mer length, while the L in the second order terms is the 
length of the actually observed filament in the tube. In 
a monodisperse solution as in our theory these quantities 
are identical. 

The importance of the second order term depends heav- 
ily on the nature of the polymers making up the network. 
In Fig. [5] the relative tube width correction obtained by 
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Fig. 6. Relative correction obtained by the second order term 
of the tube diameter (|38l) for different biopolymers as a func- 
tion of mesh size £. fd-viruses (L ~ 0.9/xm, l p ~ 2.2/^m 26 ) 
show a large correction due to their small length compared to 
the mesh size, while this effect is rather small in microtubules 
(L w 50^m, l p w 5000>m [271120] '). The correction for F-actin 
has been plotted for different length from a typical length dis- 
tribution. 

the second order term is displayed for several semifiexi- 
ble biopolymers as a function of mesh size £. It is inter- 
esting to note that the intuitive dependence on the rel- 
ative persistence length l p /L present in the second order 
term of the deflection length is rather negligible. The most 
dominant effect of the correction term is not obtained for 
the stiffest biopolymer, a microtubule, but for the small 
fd-virus. This is due to its small length to mesh size ra- 
tio. Finite length effects will influence a large fraction of 
the polymer strand and not only the boundaries. Given 
a proper control of polymer length, this effect should be 
experimentally observable in F-actin solutions. 

Focussing back on F-actin, Fig. [7] displays the result 
of our model in comparison to experimental data [28 , 29J . 
While theoretical and experimental results are certainly 
qualitatively comparable, a more detailed discussion is dif- 
ficult due to the large fluctuations of the measurements. 
However, it seems reasonable to interpret these measure- 
ments regardless of their ambiguity as an upper limit to 
the tube diameter. Two main reasons cause an experimen- 
tal observation of tube widths systematically higher than 
in the presented theory: from a technical point of view 
the microscope resolution broadens the observed tubes. 
Additionally, this effect is further enhanced by collective 
fluctuations of the complete elastic medium that remain 
unaccounted for in our approach. Contrary, the computer 
simulations presented below, can be tailored to avoid these 
effects and study the exact model system used by the the- 
ory. 



5 Simulations 

We have conducted intensive numerical simulations of the 
model system for several reasons: on the one hand they 
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Fig. 7. Comparison of tube diameter from theory, numerical 
simulations (squares) and reanalyzed experimental measure- 
ments (triangles) from [28II29] . While Dichtl has directly mea- 
sured potential strengths, Kas has recorded the maximal tube 
width a. Therefrom we estimated a lower boundary of a = a/6. 



anisotropics in performing ensemble averages. We there- 
fore choose to assume orthogonal intersections of the ob- 
stacle polymers with the plane of simulation. 

Having defined a suitable conversion for the polymer 
density, we will proceed to implement the other contribu- 
tion to the test polymer's Hamiltonian in our simulations, 
i.e. the persistence length. The test polymer is modelled as 
a chain of N rigid segments that approximate the continu- 
ous contour of a worm-like chain. The joint angle between 
two segments gives rise to a bending energy summed over 
all bonds: 

jV-1 

/3H{{U}) = kJ2 Utt+i , (41) 

i=l 

where the ti are the tangents and k is chosen such to repro- 
duce the energy of a semiflexible polymer of persistence 
length l p . The relation in two dimensions is computed to 



h(k) 



(42) 



serve as a tool to verify the validity of several approxima- 
tions used in the theoretical description developed above, 
being for example the harmonic description of the tube po- 
tential or the assumption of a single fluctuation amplitude 
for the obstacles. Furthermore, the comparison between 
the simulated transverse fluctuations and the final result 
of our theory can prove if we have succeeded in correctly 
predicting the tube diameter in a network of semiflexible 
polymers. Finally, the simulations give us the chance to 
analyze observables that go beyond the analytical theory 
presented. These are in particular distribution functions 
and open up a further possibility to comparison with ex- 
periments. 

We use a Monte Carlo simulation of a single polymer 
in two dimensions that is surrounded by point-like obsta- 
cles. This reduction will result in an equal fluctuation am- 
plitude as in the 3D model, because we have assumed the 
fluctuations along the different coordinates to be indepen- 
dent. Simulating a test polymer in 2D and measuring its 
transverse displacement, will thus on average correspond 
to either L± iV or Lj_ jZ given that the number of obstacle 
points has been chosen correctly. We calculate this num- 
ber as the number of stiff rods that cut an arbitrary unit 
area plane if these rods are of length L, density per unit 
volume v and equally distributed both in position and ori- 
entation. The approximation as rigid rods is justified by 
the large persistence length compared to the mesh size. 
The relation between polymer concentration v and point- 
density in simulations pmc then yields (|45| : 



Pmc = —vL 

7T 



(40) 



Of course, the obstacles will cut the plane under different 
angles. These can be incorporated via different statistics 
of the obstacle fluctuations. However, simulations show 
that no significant differences compared with orthogonal 
cuts occur. This can be explained by an averaging out of 



with Iq as modified Bessel function of first kind and its 
derivative Ji . 

The simulations start from an equilibrium conforma- 
tion of the test polymer and with obstacle centers r° that 
are uniformly distributed. At this point, obstacles are dis- 
criminated into those in the half-spaces left and right of 
the test polymer. During the following evolution of the 
system, every move of an obstacle and every conforma- 
tion change of the test polymer is rejected if it would 
result in a reclassification of any obstacle into the other 
half-space. Besides this constraint the evolution is only 
governed by the bending energy of the test polymer and 
a harmonic potential U(r.i) = a/2 (r^ — r°) 2 for every ob- 
stacle. During the evolution the transverse displacements 
L± of every bond from the average contour are recorded 
over the whole evolution. To avoid boundary phenomena, 
this is only done in the bulk. The whole procedure is then 
carried out repeatedly for different initial sets of random 
obstacles and random test polymer conformations. If the 
computation is repeated for different values of cr, a func- 
tion L± (tr) is obtained from which the point of self consis- 
tence Lj_(<j) — a and its error can be deduced graphically. 
Repeating the procedure for different parameters, holds 
results for the tube diameter in dependence of persistence 
length l p and concentration vL and can be compared to 
the theoretical prediction and the available scarce experi- 
mental data. As displayed in Fig.[7]the simulation results 
and the theoretical prediction to both first and second or- 
der agree remarkably well. On the basis of the available 
data any discrimination between first and second order 
would be bold. However, it has to be considered that any 
deviations due to lack of simulation time or shortcomings 
in the Monte Carlo moves will tend to reduce the observed 
tube width. The obtained simulation results are thus a 
lower boundary to the real tube diameter. 

Even if the good agreement between the theoretical 
predicted tube diameter and the values observed in nu- 
merical simulations suggests our theoretical description to 
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be valid, we employ the developed algorithms to explic- 
itly check on some of the assumption made in the course 
of deriving the tube diameter. 

One central assumption in the realm of the tube model 
is the substitution of an ensemble of neighboring poly- 
mers by an effective tube potential. This tube potential is 
modelled by an harmonic function of strength 7 as in the 
Hamiltonian |T]). This harmonic assumption seems sensi- 
ble and is also supported by preliminary experiments with 
colloidal probes [29] . Our numerical simulations can pro- 
vide further proof to the exact form of the potential. To 
this end, we have monitored the transverse displacement 
as a function of arc length. In the resulting histogram - 
see Fig. [8] (top) for some examples - we identify the dis- 
tributions maximum as the center position and analyze 
the form of the potential. Evidently, the resulting pro- 
files in the test polymer's bulk are reasonably Gaussian 
shaped, while deviations at the boundaries (compare data 
for s = 0.08 in Fig. [5J (top)) occur but are negligible for a 
tube model where L > L±. For a quantification the ratio 
of fourth moment to square of second moment of trans- 
verse fluctuations 



Q = 



(y - vY 
(y - y? 



(43) 



was considered. For a perfect Gaussian distribution this 
quantity evolves to Q — 3. As shown in Fig. [5J (bottom) 
this value is also asymptotically obtained in the simula- 
tions after sufficient simulation time. These results clearly 
support the validity of a harmonic tube potential. 

In contrast to the classical picture of an Edwards-tube 
with a rather homogeneous diameter the simulations re- 
veal a rather large variability in the local tube diameter 
as has also been observed experimentally [251125] . Carry- 
ing out extensive simulations in a large number of differ- 
ent obstacle environments allows one to record the distri- 
bution function of the tube diameters. This is of crucial 
importance, as our theoretical description has assumed 
that the tube diameter - and hence due to self-consistency 
also the obstacle fluctuation width - can be described by 
a single characteristic value. This approach only seems 
feasible if the distribution described by the characteris- 
tic value is reasonably well peaked. The simulations prove 
that the resulting distribution is indeed equipped with a 
well-defined peak (Fig. [9]). However, the variability of the 
observed tube diameters is rather large with a half-width 
of the size of the average tube diameter itself. We observe 
a sharp cut-off for small tube widths while the distribu- 
tion's tail to wide tubes is longer. The behavior at small 
tube width is dominated by the energy cost of confining 
a polymer into an increasingly smaller tube and can thus 
be considered as a polymer property. On the contrary the 
distribution at tube widths larger than the average diame- 
ter is due to void spaces. These will follow an exponential 
decay and are therefore a characteristic of the network 
architecture. 

Finally, the numerical simulations provide a means 
to explicitly check if the self-consistence is guaranteed in 
spite of the simplifying assumption of a single fluctuation 
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Fig. 8. (top) Distribution of transverse excursions at different 
arc-lengths s shows a Gaussian potential profile with rather 
large variability in the potential width. At the test polymer's 
boundaries deviations occur, (bottom) After sufficient simula- 
tion time the ratio Q (solid) approaches the characteristic value 
Q — 3 of a Gaussian distribution. The transverse fluctuation 
area (dashed) converges likewise. 
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Fig. 9. Distribution of L± sampled over polymer arc-length 
and different obstacle environments. Distributions are well- 
peaked and exhibit longer tails at large tube widths. The noise 
at diameters far from the distributions maximum is an artefact 
from the numerical discretization. 
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width. To this end, we have used the resulting histogram 
of tube diameters from above to compute a normalized 
distribution function. The fluctuation width of the obsta- 
cles are now initialized according to this very distribu- 
tion. The resulting histogram of tube diameters is then 
again fed back into the simulation as obstacle fluctuation 
distribution. This procedure is carried out until both dis- 
tributions converge against each other in a self-consistent 
manner. Surprisingly, this is already the case after the first 
iteration step of the process as displayed in Fig. [9] This 
gives strong evidence that due to the self-averaging over 
obstacles the modelling of a network with Gaussian tube 
profile and a single average tube diameter is sufficient to 
describe the physical reality. 



6 Conclusion 

We have presented a new approach to determine the ab- 
solute value of the tube diameter in semiflexible poly- 
mer networks supported by computer simulations. To this 
end the deflection length of a polymer in a hypotheti- 
cal harmonic tube was connected to the tube's diameter 
via the free energy cost for finite length polymers. The 
assumption of a harmonic tube was confirmed by sim- 
ulation results. By decomposition into independent stiff 
rods of appropriate length, we were able to establish an 
implicit equation for the deflection length. The resulting 
tube width Lj_ is in agreement with the established scal- 
ing law L± — c£ 6 / 5 Z p with mesh size £ = yJiJvL and 
persistence length l p . Our theory provides a prefactor of 
c w 0.32 and a higher order term that accounts for finite 
length effects and scales with £, 2 /L. 

The available experimental data is consistent with our 
predictions. However, its quality does not allow for de- 
tailed comparison. To provide a precise validation, we have 
complemented our theoretical work by extensive Monte 
Carlo simulations of a test polymer in an environment of 
obstacles. The resulting self-consistent tube widths per- 
fectly match the theoretical value predicted. This strongly 
supports the validity of the absolute value for the concen- 
tration dependent tube diameter. 

Furthermore, we have employed simulations to observe 
properties beyond the analytical theory. We have recorded 
the distribution function of tube widths in a network for 
different concentrations. Thereby we were able to explic- 
itly confirm self-consistency of the simplifying model with 
a fixed tube diameter. 

Both our theoretical predictions, e.g. the finite length 
contributions to the tube diameter, and our simulation 
data, e.g. the distribution functions, provide the oppor- 
tunity of feasible experimental comparison. On the theo- 
retical side, the significance of correlations and collective 
fluctuations of the complete medium, as well as an analyti- 
cal model of distribution functions may open up promising 
continuations of this work. 
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A Rigid Rod Statistics I 

To relate the polymer concentration of a network to the 
obstacle density per unit area in the simulation, we cal- 
culate the number of randomly distributed and oriented 
stiff rods per unit volume that intersect with a unit plane. 
Every intersecting rod will be described by its polar and 
azimuth angle relative to the unit plane, the point of inter- 
section and the distance between center of mass and inter- 
section point. Because of rotational symmetry the problem 
is independent of the azimuth angle and because of uni- 
form density it is independent of one of the coordinates of 
the intersection point. Hence, the problem is equivalent in 
two dimensions to the number of rods per unit area that 
intersect a unit line (see Fig.Qj^left)). The plane contains 
v rigid rods per unit area with random orientation a and 
center of mass position. The number of intersections pmc 
with the unit line (bold dashed) is computed by parame- 
terizing the center of mass (C) by the coordinate z of the 
intersection point (P) with the unit line, the distance s 
between C and P and the angle a. pmc is then obtained 
as the integral over all possibly intersecting rods: 



PMC 



da 



(44) 



,1/2 


,L/2 , 


/ dx 


/ dS 


1-1/2 J 


-L/2 JO 



where the factor 2 accounts for the fact that any rod con- 
figuration can be realized by two angles a since the rods 
have no direction. The integration area has to be cho- 
sen appropriately to only include intersecting rods. As 
r = (sin(a)s, x — cos(a)s) the Jacobian determinant of 
the coordinate transformation to integration variables is 
dr/d(x, s) = sin(a) and the integral evolves to 



B Rigid Rod Statistics II 

We derive the radial density of obstacles that effectively 
hinder the test rods fluctuations. To this end we consider 
the test rod to be aligned along the z-axis without loss 
of generality (see Fig. [10] (right)). As criterion for effec- 
tive obstruction of transverse fluctuations between test 
rod and obstacle, we demand that the line connecting 
their points of closest approach (PO) is orthogonal to 
both polymers. As a projection of the obstacle to the 
plane spanned by the test rod and PO± (dashed) re- 
covers the setup discussed in Appendix [XI the coordi- 
nates of the center of mass (C) can readily be extended 
to three dimensions by the radial distance R and the an- 
gle (3 to r = (sin(a)s sin(/3) — cos(/3)i?, sin(a)s cos(/3) + 
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sin(a)i?, x — cos(a)s) and with the Jacobian determinant 
\dr/d(x, s, R)\ = sin(a) the integration gives 



P (R) 



2v 
An 

7T 



fL/2 />7T /"27T 

dx I ds da dj3 sin (a) 

-L/2 JO JO 



-vL l 



(46) 



Consequently, the density seen by an stiff segment of length 
L in the IRM will be 



L2 2 



(47) 



Note that in the simplified geometry the obstacle density 
is not a complete radial density but a line density of ob- 
stacles on one of the four axes (positive and negative y 
and z axis). To recover the complete radial density, one 
has to sum over all of these. Hence the obstacle density 
on either one of the four axes has to be: 



7T vLL 



P= 2 — 



(48) 
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Fig. 10. (left) Sketch of a rod with center of mass (C) and 
orientation a intersecting a arbitrary unit line (dashed) to il- 
lustrate the calculation of the number of intersecting rods pmc- 
(right) Two rods with minimal distance R are considered to be 
mutually interacting only if their line of closest approach (OP) 
is orthogonal to both. 0± only serves to illustrate the analogy 
to the two-dimensional setup (see text). 



C Mode analysis of polymer and tube 

Using the dimensionless arc-length s = s/L, polymer con- 
formation y(s) — y(sL)/L, tube center y°(s) = y°(SL)/L 
and persistence length -q = l p /L the Hamiltonian JT]) can 
be written for one transverse coordinate as 



4/~ 



,~.0\21 



(y - y") 



(49) 



For free boundary conditions the Hamiltonian is diago- 
nalized by an orthonormal set of eigenfunctions ipk ~ 
sin(qk0 + 4>k) with m ir(k + 1/2) [5D]. Expanding both 
polymer and tube center in modes as y(a) — J2k Vki^k 
and y°{cr) = J2k Vk^k and using the orthogonality of the 



eigenfunctions allows one to write the Hamiltonian in the 
suggestive form 



Vk 



V°k 



1 + {qkhY 



4 A 
&k 



l + (qkld) 



(50) 



By comparison to (|4"9"|) one can read of the minimum of 
the confinement potential, i.e. the average tube center: 



(51) 



Additionally, this allows to write the complete trans- 
verse fluctuations as a sum over the inverse confinement 
strength of all modes: 



Li 



ds(x(s) — x(s)y 



E 



i 



(52) 



For the dimensionless function h(ld) relating the tube di- 
ameter to the deflection length follows: 



Hid) = 



1 



(53) 
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